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Abstract 

We show that the many-body wavefunction built as a permanent of localized non-orthogonal 
single-particle states can describe bosonic superfluidity. The criterium for the transition is ex- 
pressed in terms of the properties of the matrix of the overlaps between the single-particle wave- 
functions. We apply this result to study the superfluid-to-Bose glass transition in a disordered 
Bose-Hubbard model through a very simple variational wavefunction. We finally consider a fur- 
ther quantity, the bipartite entanglement entropy, which also provides a good estimator for the 
superfluid-to-Bose glass transition. 

PACS numbers: 05.30.Jp, 03.67.Mn, 64.70.Tg 
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I. INTRODUCTION 



The study of the interplay between disorder and interaction has recently received novel 
impulse by the experimental realizations of correlated and disordered systems through op- 
tical lattices. Particular interest has been devoted to trapped bosonic atoms in the 
presence of on-site disorder, which can be simulated experimentally using speckles or super- 
imposing incommensurate laser beams.js-S] The great opportunity offered by this kind of 
experiments is that not only disorder can be easily tuned but also the interaction among 
the atoms, this latter by means of the Feshbach's resonance technique. (gI These experiments 
have motivated a renewed theoretical activity on the phase diagram of the disordered Bose- 
Hubbard model UQ, an oldlieifl but still debated issue. This model is supposed to 
exhibit three different phases. Besides the Mott-insulating phase, occurring when interac- 
tion is strong and the particle density is commensurate, and the superfiuid phase, stable 
for weak^ interaction, in the presence of disorder an additional phase has been predicted to 
the so-called Bose glass. Such a disorder driven phase is supposed to be insulat- 



occur. 
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ing but, unlike the Mott insulator, compressible. In the absence of interaction, the Bose 
glass is simply the state obtained by condensing all bosons in the lowest energy eigenstate 
of the single-particle disordered Schrodinger equation, which is supposed to lie in the Lif- 
shits's tail of localized wavefunctions. This picture is however unstable to whatever weak 
interaction, since a macroscopy occupancy of a finite region in real space would be ener- 
getically unbearable. In fact, a more realistic view of a Bose glass is that of disconnected 
superfiuid droplets, where coherent inter-droplet tunneling is inhibited by the Anderson lo- 
calization phenomenon, hence the insulating behavior, yet transferring electrons between 
different droplets is costless, hence the finite compressibility. 

In this paper, we investigate the superfiuid-to-Bose glass transition by means of a very 
simple variational wavefunction. It consists of a permanent of non-orthogonal single-particle 
wavefunctions that are determined in a variational manner. Although the approach is not 
as rigorous as e.g. quantum Monte Carlo, jsl, [lo|, 13, 15| nevertheless it provides a very 



transparent description of the transition and of the same Bose glass. Indeed, all coherence 
properties of the wavefunction are hidden in the single-particle wavefunctions that build the 
permanent. In particular, since these wavefunctions need not to be orthogonal, it is their 
overlap matrix that seems to play an important role. Actually one can derive an approximate 
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criterium for the permanent to have macroscopic condensation at zero momentum, hence 
superfluidity, which involves just the overlap matrix. We find that this criterium agrees well 
with the more rigorous one based on the superfiuid stiffness. 

Finally, in the last part of the paper we study the single site von Neumann entropy, 
which could also be used to identify the superfluid-to-Bose glass transition. The quantum 
entanglement applied to many body systems has attracted lot of theoretical interest in 
recent years (see for instance Refs. 22h25| and references therein). So far, however, the 
entanglement witnesses have never been used in the context of disordered Hubbard model. 
We show that, within our variational approach, the Bose glass seems to be characterized by a 
finite probability of the single-site von Neumann entropy to have zero value. In other words, 
as disorder increases, the entropy distribution gets broaden and broaden until the probability 
to have zero entropy becomes finite, which we identify as the onset of Bose glassiness. This 
criterium is qualitatively good as far as we deal with small system size since the zero entropy 
is the signature of vanishing charge fiuctuations in some regions of the sample which becomes 
then easily disconnected. A more quantitative analysis for large systems would involve the 
study of the entropy distribution and the appearance of the bimodal profile in agreement 
with the droplet scenario. Our conjecture is that the weights and the positions of the two 
modes at the transition should be related to the percolation threshold. 

The paper is organized as follows. In Sec. |TT]we present a criterium for the onset of super- 
fluidity that is based on the matrix of the overlaps between the single particle wavefunctions 
that built the variational many-body wavefunction. In Sec. IIIII we introduce the disordered 
Bose- Hubbard model and in Sec. [IV] we describe the variational many-body wavefunction 
that we study. Sec. |V]is devoted to studying the bipartite entanglement entropy as another 
tool to identify the superfluid-to-Bose glass transition. Final conclusions are summarized in 
SeclVD 



II. OVERLAP MATRIX CRITERIUM 



Let us imagine a Hartree-Fock-like variatonal approach to the problem of disordered 
and interacting bosons that amounts to search for the permament that minimizes the total 
energy. Unlike conventional Hartree-Fock theory for fermions, such an approximation for 
bosons does not lead to signiflcant simplifications with respect to exact numerical simulations 
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because permanents are extremely difficult to handle with. Nonetheless, let us assume we 
have successfully performed the calculation and found the optimal permanent for particles, 
which can be written as 

a=l \ i J 

in terms of a set of single particle wavefunctions ipa, with amplitude ipai at site i. In ([T]) 
bl creates a boson at site i. Unlike Slater determinants, the wavefunctions that built a 
permanent need not to be orthogonal one to each other, so we expect the overlap matrix Q 
to have non-zero off-diagonal elements flai3 = '^ai^*pi- instance, if all bosons condense 
into a single state, then all ilja are equal and VLap = 1, Wa,f3 G [1,A^]- In the generic case 
where ipa are distinct, we may wonder whether wavefunction ([T]) describes a superfiuid. In 
what follows we derive a simple criterium that is based on the properties of the overlap 
matrix. 

One can verify that the norm of ([1]) can be written as an integral over classical variables 



as 
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(^|v^)= / JJ ^M^e-^(«'«^) (2) 



with the following action (see Appendix |A]) 

a 

We assume that the main contribution to the integral comes from the saddle point, i.e. the 
solution of 

= E ^ 1 , |2 ^P- (4) 

The above equation implies that finite values of C,a appear in groups, or equivalently that Q 
is a block matrix. In the presence of disorder this is suggestive of the existence of clusters 
occupied by bosons whose wavefunctions mutually overlap. Because of interaction, a cluster 
can not accomodate all particles unless it covers all the system, which would correspond also 
to superfluidity. If we linearize ([1]), we find that the condition for the apperance of a cluster 
reads 

which corresponds to a block of Q that acquires an eigenvalue greater than two. As we 
mentioned, this is still not the condition for superfiuidity. The latter rather implies that a 



block in fl should grow, or several blocks should merge, i.e. start to overlap, till a percolating 
cluster emerges. This condition is likely to be equivalent to an eigevalue of Q that grows 
with the number of bosons A^. We finally mention that the saddle point approximation that 
we used is rigrously valid only if blocks are big enough. 

A. Density matrix and overlap condition 

A better and more transparent criterium to detect a long range order is to resort to the 
definition of the density matrix. 



Off-diagonal long-range order implies that Cij is finite for |« — j| — > c>o. Within the path 
integral formulation, using Eqs. (lA2p . (IA6p . one can verify that Cij can be written as 



a/3 

where (...) is now the average weighted by e""^, with the action S given by Eq. ([3]). If the 
saddle point of the action is characterized by finite ^a, one could be tempted to set 



This is not fully correct. Indeed, if f2 is a block matrix, within each block only the relative 
phases of the C,a are fixed, while an overall phase is still free and has to be integrated out. 
This implies that is correct only if a and /3 are within the same block, otherwise the 
relative phase between the two blocks will suppress the average. This suggests that off- 
diagonal long-range order sets up only if a single block percolates. More rigorously, we 
define 



a, = {^\blb^\^). 



(6) 




(7) 



(8) 




(9) 



a 



through which the density matrix reads 



(10) 



The saddle point equation in terms of is 




(11) 
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where 



(12) 



a 



is the density matrix of the wavef unctions. An extreme superfluid solution identified by 
^i — ^, Vi could be stabihzed if 



where Ng is the number of sites. 

B. An example: the bosonic crystal 

As a simple application of the previous results we shall now investigate the possibility 
that a permanent that describes at the "Hartree-Fock" level a Bose-Wigner crystal, could 
be also superfluid, actually a supersolid. Let us consider a commensurate density N/Ng < 1 
of interacting bosons on a lattice with A^^ sites in the absence of disorder. If the repulsion 
is sufficiently strong and long ranged, we may imagine that the best variational permament 
wavefunction describes a bosonic superlattice with Bravais vectors 



where a is the superlattice parameter. We write the Wannier single-particle wavefuncions 
that correspond to the Bose-Wigner crystal as 



where k runs within the Brillouin zone of the original lattice, Rj spans all lattice sites while 
R only the superlattice ones. The permanent is therefore 




(13) 



R = a(mi,m2, ...,md), 



(14) 




(15) 




(16) 



R 



In this case the overlap matrix is 




(17) 



We assume for simplicity that the wavefunction is gaussian. 




(18) 
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where v is the volume of unit cell and i the localization length, so that the overlap matrix 
has the simple expression 

2 |R-R'|2 

fiRR, = e-"T, (19) 

and depends only on the distance, Qrr> = ^^r-r'. From the saddle point equation we find 
that the condition to have superfluidity is simply 

^^]R = ^e-"'^ >2. (20) 

R R 

By means of the Jacobi theta function 6'3(0|x) = ^^=-00^™^ ^'i- dSO]) means 

03 (0\e--"-"/'"^ > \/2. (21) 

The condition Eq. fl2T|) fixes the critical overlap between the Wannier functions, parametrized 
by the ratio i/a, above which bosons condense at zero momentum, i.e. the many-body 
wavefunction describes a supersolid. 

In the next section we shall consider the Bose-Hubbard model in the presence of disorder 
which causes now the localization of the single-particle wavefunctions. We are going to see 
that, also in that case, in spite of the Anderson localized nature of the single-particle states 
obtained variationally, the many-body wavefunction can be superfiuid. 



III. THE MODEL 



We consider a system of interacting bosons on a disordered 2-dimensional lattice with 
Ns = L"^ sites, described by the following Bose-Hubbard hamiltonian 



rii^Ui - 1), 



(22) 



where (b-) creates (annihilates) a boson at site R,, {ij) denotes the sum over all pairs 
of neighboring sites, Ui = b\b^ is the boson local density, and, finally, are random on-site 
energies uniformly distributed between —A and A. 



16 



171] and Fisher et al. 18|, the Hamil 



Since the seminal works of Giamarchi and Schulz 
tonian (l22l) has been studied with a variety of techniques, mainly in one and two dimensions. 
More recently, highly sophisticated numerical simulations {tI. IsI. Iiol. 
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2l| have been 



performed to uncover the full phase diagram and settle up some debated issues, like the 



possibility of a direct superfluid-Mott insulator transition. The results we are going to 
present are by no means comparable in accuracy with the aforementioned numerical sim- 
ulations. Our main purpose is not to compete with those simulations, but just to provide 
an interpretation of the phase diagram in terms of a simple Hartree-Fock-like single-particle 
picture. 



IV. THE METHOD 



The simplest way to deal with interacting electrons is the Hartree-Fock approximation, 
which amounts to search for the best wavefunction within the subspace of Slater determi- 
nants. This approximation reduces the complicated many-body problem to a single-particle 
one with an effective potential generated by all other particles. Even more realistic ap- 
proaches, like the Density- Functional theory in the Local Density approximation, eventually 
reduces to the self-consistent solution of a single-particle Schrodinger equation. The great 
advantage is that a single-particle description is very intuitive and, although it could be too 
naive in many cases, at least it is a simple starting point for more complicated approaches. 

The obvious generalization of the Hartree-Fock approximation to interacting bosons 
would be searching for the best wavefunction within the subspace of permanents, the bosonic 
analogues of Slater determinants. However, as we already mentioned, unlike Slater deter- 
minants, permanents are well defined even if the single-particle wavefunctions that are used 
are not orthogonal to each other. Therefore the optimization procedure is not reduced to 
solving a single eigenvalue problem, which would produce a set of orthogonal wavefunctions, 
but becomes rather complicated, practically unfeasible, hence further approximations are 



required. 
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Our simplified approach consists in adding one boson at a time with a wavefunction that 
is the ground state of a non-interacting Hamiltonian with a potential determined by the 
already added bosons. Specifically, we consider the not normalized N boson wavefunction 

\'^N) = f[[J2Xa,bn\0) (23) 

Q=l \ j / 

where Xaj, j = l,---,-?^^ are generically non-orthogonal single-particle wavefunctions. The 
first wavefunction xij is the ground state of (l22l) with U = (no bosons are present). The 
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(M + l)-th wavefunction is instead the ground state of 

t 



Happ — 2 



J2^lb,+J2['^ + U{n.^'^n,, (24) 



where 

is the average density of the previously added M bosons. We define an M x M matrix 

Ns 

= ^o.0i% (26) 

i=l 

a, /3 = 1, . . . , M, and, for each couple of lattice sites (i, j), the (M + 1) x (M + 1) matrices 



where Xi = (xij? •••) XmiY- It follows that the mean local density that is required for adding 
the next (M + l)-th boson, namely Eq.( l25|) . can be written as follows 

<->-^- (-) 

where Per(X) is the permanent of X. At each iteration, one can also calculate the inter-site 
density matrix, since 

_{^gb^_FeriD^ 

= " ~ ^''^ ^^^^ 

hence investigate the eventual offset of long-range order. This procedure is iterated until the 

desired number of bosons is reached. We note that the method can be easily extended to 

study excited states - it is sufficient to select at any iteration not the ground state but an 

excited one - hence at finite temperature, even though in what follows we just focus on the 

lowest energy states. 

The superfluid properties of the model can be accessed by calculating the superfluid 
stiffness psf defined through 

r:^ fi'^ R„ 

(29) 



L2 d^Ee 

Psf ' 



N 

where Ee is the average value of the Hamiltonian ( !22l) with twisted boundary conditions 
along a given direction x, or, alternatively, with hopping parameter between neighboring 
sites tij = te*^*"'^, where r = (R^ — Rj)/L. In terms of permanents we calculate 



-.^-E(|-M..--))c..4E(?if-) 



(30) 



where Cij is defined in Eq. ( 128|) and Veii^Iijki) 

I 



with 



(31) 



^ Xi Xj 

'-ijkl = xl Sik Sjk 

\xl Sil Sji J 

Upon repeating this calculation for several disorder configurations, averaging over them and 
setting equal to zero the stiffness in the region of parameters where its variance is greater 
than its average (cutting, therefore, the values of psj statistically undetermined), we finally 
obtain the phase diagram of Fig{T|, which is a contour-plot of the averaged superfluid stiffness 
Psf of a 2-dimensional model with filling fraction u = N/Ng = 1/4. We note that in a finite 
system there cannot be a true gauge-symmetry breaking, hence the phase diagram Fig{T] 
is just an indication of what could happen in the thermodynamic limit. Nevertheless, the 
qualitative behavior that we find is physically sensible: psf decreases on increasing disorder 
and, at fixed disorder, first increases with U and then diminishes. In Fig. [T] and in the 
following figures the values of U and A are in units of t. 
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FIG. 1: (Color online) Contour-plot of the superfluid stiffness p^f, Eq. (j29p . for a 6 x 6 square 
lattice at filling = 1/4, averaging over 400 disorder configurations. The two thick dashed lines 
and the dotted line show the border of the superfluid phase, by looking at some indicators related 
to the overlap matrix, as discussed in Sec. IIV A[ 



We previously mentioned that a realistic view of a Bose glass is that of disconnected 
droplets. A way to confirm this idea is plotting the average density as shown in Fig. |2l 

10 



where we used, instead of Eq. fl2^ . (rij) ~ IXaiP? which is a good approximation for 
local densities, in order to computationally reach larger system sizes. We find that, at 
large disorder, bosons are indeed concentrated into droplets whose magnitude increases on 
decreasing disorder until a percolating cluster appears, which must presumably signal the 
onset of superfluid. 




FIG. 2: (Color online) Contour plot of the particle density (rij) on a 28 x 28 lattice with = 49 
bosons, i.e. filling v = 1/16, at [/ = 10 and for a single realization of disorder with A = 1,2,3 
from right to left (i.e. crossing the transition, see inset of Fig. E]). The darker the sites the higher 
is the density. 

A. Comparison with the overlap matrix method 

Let us compare the transition obtained by the superfluid stiffness with that given by 
the overlap matrix criterium. In our case, given the many-body wavefunction Eq. (!23|) . the 
overlap matrix VL is that defined in Eq. fl26|) with ipai = Xai- The action is given by Eq. ([3]) 
and the saddle point equation reads as in Eq. ([5]), which has non trivial solution if Q has 
eigenvalue 2. 

In Fig. [T]we plot a dotted line below which {Q — 21), or equivalently, (0 — 21), has both 
positive and negative eigenvalues for any disorder configuration, implying that the saddle 
point equation has always non trivial solutions. Above that line, instead, for some disorder 
configurations all the eigenvalues of Q are smaller than 2. 

On the same figure. Fig. [H we plot two thick dashed lines which correspond to J-" = 2 
(the upper line), namely, when the value of J-" as defined in Eq. (fT3ll . averaged over 400 
disorder configurations, equals 2, and min{J^) = 2 (lower line), namely, when the minimum 
value of J-" among the disorder configurations equals 2. Below those lines the corresponding 
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quantities exceed the threshold value. Notice that a non trivial solution of the saddle point 
equation occurs also for large disorder and weak interaction (the dotted line keeps on growing 
decreasing U) while J-", which detects the long range order (dashed lines), follows correctly 
the stiffness behavior also for small interaction. 

In summary, we find that the superfiuid properties of the many-body wavefunction f l23p 
can be related to the overlap matrix Q between the single-particle wavefunctions Xai- This 
result also clarifies why, even though each wavefunction Xai is the lowest energy solution 
of the Schrodinger equation of a particle in a disordered potential, hence would be always 
localized in two dimensions, and also in higher dimensions if, as presumably is the case, it 
lies in the Lifshitz tail, nevertheless the permanent built with them could still be superfiuid. 

A final comment. We optimized the wavefunction by explicitly evaluating permaments. 
This is in general very cumbersome, not much simpler than an exact numerical solution of 
the problem, which is the reason why our simulation size is small. However, we could adopt 
an oversimplified approach and evaluate the Hartree potential in Eq. (l24l) for the (M+ l)-th 
boson using 

M 

(n,) ^ l^-l'' (32) 

a=l 

as if the already present M bosons were distinguishable. This approximation simplifies a 
lot the procedure to determine the single-particle wavefunctions, which can be pushed to 
very large system sizes. These single-particle wavefunctions can then be used to construct 
the permanent, whose superfiuid properties can be assesed by the overlap matrix criterium. 
Even though such a procedure is hardly justifiable from the variational point of view, it 
provides a phase diagram qualitatively correct. Using this approximation one can study the 
spectrum of the overlap matrix Q for larger system, finding that, indeed, the transition is 
characterized by the fact that the greatest eigenvalue of Q starts growing towards value A^, 
the number of bosons, for A going to zero. From Fig. [3] one can see clearly, in fact, that a 
single eigenvalue of Q (although quite noisy, it is a single level) separates from the others at 
A 4, approximately at the same value of disorder as that obtained for a smaller system size 
with (rij) given by fl271) (cfr. the transition line at f/ = 10 in Fig. [1]). The other eigenvalues 
(see the inset of Fig. [3]) accumulate around zero below the transition, and around one above 
it, at least for filling not greater than one, as in our case. 
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FIG. 3: Eigenvalues of ft as functions of A, at a given configuration of disorder, for 49 bosons on 
a square lattice with 256 sites, namely always for filling 1/4, at U = 10. The inset magnifies the 
lower part of the main plot. 

V. SPATIAL ENTANGLEMENT ENTROPY 

Another quantity which may be interesting to look at is the single site entanglement 
entropy. We define as Pnii) the probability to have n bosons at site i, which must trivially 
satisfy '^^=0 Pn{^) — 1) where N is the total number of bosons. The single-site entropy Si 
is thus defined through 

N 

Si = - ^ Pn{i) lnp„(i). (33) 

n=0 

In a disordered system it is also convenient to define its probability distribution through 

^('^) ^disorder, (34) 

which is obtained considering all sites and all disorder configurations. In spite of its simple 
definition, we are going to show that the single-site entanglement entropy, and especially its 
distribution probability, is a quite enlighting quantity in the presence of both disorder and 
interaction. 

A. A limiting case 

We start considering the case where all bosons condense into a single state, i.e. Xaj — Cj^ 
Va, with the normahzation condition '^fliQ — 1- In this case Per(fi) = N\ and the state 
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Eq. fl23|) . including the normalization factor 1/ y^FeiiJV) can be decomposed in the site Fock 
basis as 



i^) = ;7= EcaM io) = Ew^nC(^lr^|o)- (35) 

V i / {n} * i=i 

The sum runs over {n} = (rii, ?t,2, "".atJ, the configurations of occupation numbers with 
the constraint Yli = We calculate the reduced density matrix by partitioning the 
sites in two blocks: [1, i] and [i, Ng], and tracing out the sites belonging to the second block: 

p'= E — ^(oi n n (36) 

n<.+i,..,njv3 lli=£+l '^^^ i=i+l i=l+l 

This is a diagonal ^i^^j^^/l^^; x pj^^^l^^fy^ matrix. For each diagonal element Pni,...,ne, corre- 
sponding to the configuration of the occupation numbers (ni,n2, ...yni) of the i sites that 
are not traced out, we obtain the following expression 



From now on we shall focus on the simplest partitioning, keeping only one site and tracing 
over all the others, i.e. i = 1. In this case we get a very simple binomial expression of the 
reduced density matrix p„ = Pn(l) 

p.= (^^j (i-iciH^'^-^^icir (38) 

It is straightforward to check that XlnP" ~ 1- The reduced density matrix is fully local; 
it depends only on the value of the wavefunction ( on that site. We can now calculate the 
entanglement von Neumann entropy at site 1, 5*1, through fl5^ . From Eqs. (1381) and fl55]) 
we find that for any fixed value of |Ci| G (0, 1) and for large A^, the asymptotic value of 5*1 is 

^1 > Sioc = l\nN + A, (39) 

with A = I {l + In [27r|CiP(l - KiH] }• At the two extremes, |Ci| = 0, 1, we have Si = 0. 

Let us consider first the non-interacting case with disorder. All bosons condense into a 
localized wavefunction. Within the localization region Si ~ In A^, while outside Si = 0. We 
conclude that the probability distribution ( l34l) becomes peaked at the single value 5 = 
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in the thermodynamic hmit, Ng ^ oo, where the locahzation region has zero measure with 
respect to the whole system. 

On the contrary, without disorder and deep in the superfluid phase, small U, we expect 
that 

lOP^^^^ (40) 

where the filling fraction u = N/Ng. In this case, the single-site entropy in the thermody- 
namic limit is finite, site-independent and depends exclusively on the filling fraction u: 



oo 



Si ^ 5sF = -Ini^) + 6"'^ V— Inn! (41) 

N^oo ^-^ nl 

n=0 

The sum Eq. pil) converges for any value of u. As one can see from Fig. |U the thermody- 
namic limit Eq. (14T]) is approached already with few bosons, for small u. In the same figure 
we also draw the maximum entropy line, Eq. fl39|) . which indeed lies above all curves. 

Si 
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FIG. 4: (Color online) Single site entanglement entropies obtained by Eqs. (j38p . (j33p . Upper 
dot-line: the maximum entropy, reached when = 1/2. Upper middle dot-line: entropy for 
the superfluid with v = 7/2 case, i.e. for \C,i\^ = 7/{2N). Lower middle dot-line: entropy for the 
superfluid in the commensurate 1^ = 1 case, i.e. for = 1/A^. Lower dot-line: the superfluid 
entropy for |CiP = l/(4iV), i.e. u = 1/4:. 

In conclusion, we expect that the probability distribution of the single-site entropy is 
peaked at a finite value deep in the superfiuid phase while at value zero deep in the Bose- 
glass phase. The obvious question is what happens in between. 
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B. Single-site entropy across the superfluid-to-Bose glass transition 



Let us now consider the state Eq. (123!) . including the normahzation factor 1/ A/Per (O), 



which in the Fock basis of sites reads 



1 



VI/) 



(42) 



M {fe}}i=i"=i 



The sum runs over {n} = (ni, n2, njvj, the configurations of occupation numbers with 



the constraint Ei^* = iV, and over fe} with = Efr^} Efei^^} Etel^te},^} - 



where Pj{a) is the particle label which takes rtj integer values among values. The reduced 
density matrix for a single site is found to be 



The explicit expressions of the reduced density matrix for N = 2 and = 4 are given in 
Appendix [Bl In what follows we will consider the simple case of A^ = 4 bosons. In Fig. [5] 
we show the probability distribution of the single site entanglement entropy , which, in the 
absence of disorder, is peaked around S ~ 0.62, broadens under the action of disorder and 
finally, when disorder becomes strong, develops again a peak but at S* = 0. We observe that 
the probability of having zero entropy becomes non-zero only above a disorder threshold 
that occurs approximately when also the stiffness vanishes, as one can see from Fig. |6l 

We have found that in our finite size simulation the crossover between the two limiting 
behaviors, P{S) peaked either at S 0, deep in the superfluid, or close to S* = 0, deep in 
the Bose glass, is continuous. We have also noticed that the point at which P{S) becomes 
finite at S* ^ seems to coincide with the point where the stiffness becomes very small, 
which could be taken as the finite size signal of the superfluid-to-Bose glass transition. This 
suggests that the behavior of P{S = 0) could be used as well to identify the superfiuid- 
to-Bose glass transition in alternative to the superfluid stiffness in any finite-size numerical 
simulations. Moreover, we note, looking at Fig. [5l that, close to the transition (A ^ 4), the 
profile of the distribution (the green curve) is almost an equally weighted bimodal, peaked 
either at small S (insulating regions) and at 




2 



(43) 



S ^ SsF{l^eff) 



16 




0.2 0.4 0.6 0.8 1 1.2 
S 



FIG. 5: (Color online) Probability distribution of entropy for f7 = 10 and for different values of 
disorder strength (A = 1,2,4,6). The plot has been made after lO'* realizations of disorder in a 
system of = 4 bosons on a 4 x 4 square lattice. In the inset: the fluctuations of S increase with 
the strength of disorder and saturate to the maximum value of entropy for 4 bosons which is ~ 1.4 

(superfluid regions) given by Eq. fHTj) where u is replaced by z/g// = ^/Pc- The numerical 
result for the second peak close to the transition (5* ~ 0.8, for u = 1/4), is consistent with 
Pc ~ 0.59 which is the site percolation threshold for a square lattice. We have checked 
this result also with other filling fractions. This finding is in a nice agreement with the 
percolating droplets scenario. 

VI. CONCLUSIONS 

We have considered a two-dimensional Bose- Hubbard Hamiltonian in the presence of dis- 
order, and constructed a trial many-body wavefunction by solving a single-particle problem 
for each boson at a time that feels the effect of all the others as an effective potential. 
By this wavefunction we have studied the transition between the superfluid and the Bose 
glass. Beside the vanishing stiffness, we find that the transition can be characterized by the 
eigenvalues of the matrix of the overlaps between the single-particle wavef unctions. Another 
quantity we propose to look at is the single-site entanglement entropy. In particular, we 
argue that when the probability to measure zero entropy becomes finite, P{S = 0) > 0, the 
superfluid start to vanish and the Bose glass phase sets in. The bimodal entropy distribution 
is, then, the signature of lack of percolation among superfluid clusters. 
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FIG. 6: (Color online) The minimum site entanglement entropy Smin (light-green lines) compered 
with the superfluid stiffness psf (dark-red lines) as functions of the disorder, for U = 10 and for 
filling fraction z/ = 1/4 (main plot) and v = 1/16 (inset). The stiffness, psj, is taken as an average 
among 400 configurations of disorder for a system of = 9 bosons on a square lattice of 6 x 6 
(z^ = 1/4, main plot) and 12 x 12 = 1/16, inset). The minimum entropy Smin is obtained taking 
the minimum values of S among 10^ configurations of disorder, and for A = 4ona4x4(i/ = 1/4, 
main plot) and 8 x 8 {u = 1/16, inset) lattices. The blue arrows point at the transitions, obtained 
from the overlap matrix method as explained in Sec. [Ill 

Appendix A: Derivation of the action 

Let us consider the following wavefunction 



which generalizes Eq. ([T]) where the integers are all equal to 1. Let us define also the 
overlap matrix Qai3 = Yli '^ai'^P^i which is Hermitian and positive defined, such that can 
be parametrized as = ipip'^ = V^X\'^V, in terms of a unitary matrix V. Let us suppose 
now that another unitary matrix U exists so that we can write ip = V'^XU, and define the 
following bosonic operators 




(Al) 




(A2) 




(A3) 



a 



Defining also 



(A4) 
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one can then verify that 

J2^^4 = J2vLXabi = e^bie 

i a 

where we have used the Hausdorff relation 



Xahi = e^bie-^. 



Therefore Eq. (lAip can be written as follows 

i^)=e^nw)"io) 

a 

Moreover, one can check that 



-r 



(A5) 



(A6) 



(AT) 



with |Aap = 1 + |faP- Now, in order to rewrite Eq. (lASP on the basis of ba instead of ba, we 
have to define Vai = Xla KJa^a^ai and ^a = J2ia ^ai^LK^ and noticing that vv'' = Q - 1, we 
finally have 

a /3 7 

Expanding the exponents and using the commutation relations, in particular the equality 

^m^t"-|o) = _^l_5t("--™-)|o) for m <n, one can verify that 



{'m\y{na — m)\ 



1 2m 



(AlO) 



/3 ■y a m=0 

equal to Ha using the definition of the generalized Laguerre polynomials, 

Ltiz) = yi"_n ui f""\f?' — which can be written also in terms of confluent hyper- 
geometric functions, Ij^{z) = -^^j^ n, k + 1; z). As a final result the norm of (lAip can 



be written as an integral, Eq. ([2]), with the action 



n, 



(All) 



If = 1, Va, since L\{—\^a\'^) = (1 + l^^p), the action reduces to Eq. The saddle point 
equation reads then 



/3,n^>0 



(A12) 
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which reduces to Eq. (jl]), for all = 1. Since, for <^ 1, we can expand (- 
^"1^,-" ^1 + ^^I^^P^ , the eigenvalue equation, when discarding 0{^^) terms in the r.h.s. of 
Eq. (lA12p . becomes 

which generalizes Eq. ([5]) for arbitrary n^. 

Appendix B: Reduced density matrix: two simple cases 

Here we are going to show within a toy model made of 2 bosons arranged on 3 sites how 
disorder and interaction can conspire to enhance the bipartite entanglement entropy. With 
two bosons N = 2 and without disorder, in the delocalized phase, = 1/3, from Eqs. 
fl38|) . fl33|) . we get the following single site entropy 

16 

5 = 21n3 - y ln2 ^ 0.96 (Bl) 

What we shall be seeing in the following is that this value can be overcome introducing a 
suitable amount of disorder. 

If we now take Eq.f l42p as the many body wave function, we get, for N = 2 bosons and 
generic Ns sites, the following three diagonal elements of the reduced density matrix 

Po = Yl (^2|xi^X2i|' + l^i^^2i + X2^Xlj\^j (B2) 



i=2 

2 



P2 = ^^\XnX2,\ m 

Let us now consider Ng = 3 and a simple on site disorder such that = ±A. We can have 
therefore 8 possible configurations of disorder. For [/ ^ A, by semiclassical considerations 
and using Eqs. flB2tiB4p . we can have probability ~ 1/4 to have 5* ~ 0, i.e. P{S ^ 0) ^ 1/4 
while P{S ~ ln2) ~ 1/2 and P{S ^ 0.96) ^ 1/4, the same entropy as without disorder. For 
t/ < A, instead, we have P{S ~ 0) ~ 1/2, P{S ~ 0.96) ~ 1/4 and P{S ~ | ln2 ~ 1.04) ~ 
1/4, namely we can have a sizable probability to have maximum entropy for two bosons, so 
to exceed the value in Eq. (]B1|) . The disorder, therefore, peaks the probability distribution 
at zero while widening the entropy fluctuations. 
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In general for filling u < 1, the superfluid single site entropy is S* ~ z/(l — Inu), as shown in 
the text. By introducing strong enough disorder there is the possibility for a fraction k of 
the total sites to have very large local energies which make those sites inaccessible for the 
particles. This induces a larger effective filling, j-^, and consequently enhances the entropy. 

For the sake of completeness we hereafter report the form of the reduced density matrix 
for = 4 bosons, used in the paper in Sec. IV B[ To simplify notation we define the matrix 

f xii xij xik xii ^ 

X2i X2j X2k X2I 
X3i X3j XSk X3l 

\ Xii Xij Xik Xii J 

so that the single site reduced density matrix can be written in the following form: 



£ijkl 



(B5) 



Po 



Pi 



P2 



P3 



pi 



■'iijk\ 1 2 



Ns ( Ns Ns Ns , Ns Ns 

^ ' 1=2 K j>i k>j l>k j^i k^i,j 

j>i ' jj^i 

Ns f Ns Ns ^ Ns ^ 

^ ^ 1=2 L j>i k>j ' j^i ' J 

k^Eip-(^^"^)i^ 



(B8) 



Per(^]) 3! 



Per(l]) 4! 



i=2 



|Per(£ 



1111m2 



(B9) 



(BIO) 



From equations above one can easily guess the form of the reduced density matrix for a 
generic value of N. 
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